// 版权所有2009 Go作者。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// Exp返回e**x，x的基e指数。
// 
// 特殊情况是：
// Exp（+Inf）=+Inf 
// Exp（NaN）=NaN 
// 非常大的值溢出到0或+Inf。
// 非常小的值溢出到1。
func Exp(x float64) float64 {
	if haveArchExp {
		return archExp(x)
	}
	return exp(x)
}

// 下面的原始C代码、长注释和常量
// 来自FreeBSD的/usr/src/lib/msun/src/e_exp.C 
// 并随附此通知。go代码是原始C.
// 
// ================================================================================================
// 版权所有（C）2004由Sun Microsystems，Inc.保留所有权利。
// 
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// =====================================================================
// /
// /
// exp（x）
// 返回x的指数。
// 
// 方法
// 1。参数缩减：
// 将x缩减为r，使| r |<=0.5*ln2~0.34658。给定x，求r和整数k，使得x=k*ln2+r，| r |<=0.5*ln2。
// 
// 此处r将表示为r=hi-lo，以获得更好的
// 精度。中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网中国英语学习网。用特殊有理函数逼近exp（r）在
// 区间[0,0.34658]：
// Write 
// r（r**2）=r*（exp（r）+1）/（exp（r）-1）=2+r*r/6-r**4/360+。。。
// 我们在[0,0.34658]上使用一个特殊的Remez算法生成一个近似R的5次多项式
// 该多项式近似的最大误差
// 为2**-59。在
// 换句话说，
// R（z）~2.0+P1*z+P2*z**2+P3*z**3+P4*z**4+P5*z**5 
// （其中z=R*R，P1到P5的值如下所示）
// /和
// /| 5 |-59 
// /| 2.0+P1*z+…+P5*z-R（z）<=2 
// /| 
// /exp（R）的计算因此变成了
// 2*r 
// exp（r）=1+----
// r-r 
// r*R1（r）
// =1+r+-------------（为了更精确）
// 2-R1（r）
// 其中
// 2 4 10 
// R1（r）=r-（P1*r+P2*r+…+P5*r）。中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译：中国英语学习网编译。缩小比例以获得exp（x）：
// 从步骤1开始，我们有
// exp（x）=2**k*exp（r）
// 
// 特殊情况：
// exp（INF）是INF，exp（NaN）是NaN；
// exp（-INF）是0，而
// 对于有限参数，只有exp（0）=1是精确的。
// 
// 精度：
// 根据误差分析，误差始终小于
// 1 ulp（最后一位的单位）。
// 
// Misc。信息。
// 对于IEEE双
// 如果x>7.09782712893383973096e+02，则exp（x）溢出
// 如果x<-7.4513321101941108420e+02，则exp（x）底流
// 
// 常数：
// 十六进制值是以下
// 常数的预期值。可以使用十进制值，前提是
// 编译器将十进制值转换为二进制值的精度足以产生所示的十六进制值。

func exp(x float64) float64 {
	const (
		Ln2Hi = 6.93147180369123816490e-01
		Ln2Lo = 1.90821492927058770002e-10
		Log2e = 1.44269504088896338700e+00

		Overflow  = 7.09782712893383973096e+02
		Underflow = -7.45133219101941108420e+02
		NearZero  = 1.0 / (1 << 28) // 2**-28 
	)

	// 特殊情况
	switch {
	case IsNaN(x) || IsInf(x, 1):
		return x
	case IsInf(x, -1):
		return 0
	case x > Overflow:
		return Inf(1)
	case x < Underflow:
		return 0
	case -NearZero < x && x < NearZero:
		return 1 + x
	}

	// 减少；计算为r=高-低，以获得更高的精度。
	var k int
	switch {
	case x < 0:
		k = int(Log2e*x - 0.5)
	case x > 0:
		k = int(Log2e*x + 0.5)
	}
	hi := x - float64(k)*Ln2Hi
	lo := float64(k) * Ln2Lo

	// compute 
	return expmulti(hi, lo, k)
}

// Exp2返回2**x，x的基2指数。
// 
// 特殊情况与Exp.
func Exp2(x float64) float64 {
	if haveArchExp2 {
		return archExp2(x)
	}
	return exp2(x)
}

func exp2(x float64) float64 {
	const (
		Ln2Hi = 6.93147180369123816490e-01
		Ln2Lo = 1.90821492927058770002e-10

		Overflow  = 1.0239999999999999e+03
		Underflow = -1.0740e+03
	)

	// 特殊情况
	switch {
	case IsNaN(x) || IsInf(x, 1):
		return x
	case IsInf(x, -1):
		return 0
	case x > Overflow:
		return Inf(1)
	case x < Underflow:
		return 0
	}

	// 参数约简；x=r×lg（e）+k带| r |≤ ln（2）/2。
	// 计算为r=hi-lo以获得额外精度。
	var k int
	switch {
	case x > 0:
		k = int(x + 0.5)
	case x < 0:
		k = int(x - 0.5)
	}
	t := x - float64(k)
	hi := t * Ln2Hi
	lo := -t * Ln2Lo

	// compute 
	return expmulti(hi, lo, k)
}

// exp1返回e**r×2**k，其中r=hi-lo和| r |≤ ln（2）/2。
func expmulti(hi, lo float64, k int) float64 {
	const (
		P1 = 1.66666666666666657415e-01  /* 0x3FC55555; 0x55555555 */
		P2 = -2.77777777770155933842e-03 /* 0xBF66C16C; 0x16BEBD93 */
		P3 = 6.61375632143793436117e-05  /* 0x3F11566A; 0xAF25DE2C */
		P4 = -1.65339022054652515390e-06 /* 0xBEBBBD41; 0xC5D26BF1 */
		P5 = 4.13813679705723846039e-08  /* 0x3E663769; 0x72BEA4D0 */
	)

	r := hi - lo
	t := r * r
	c := r - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))))
	y := 1 - ((lo - (r*c)/(2-c)) - hi)
	// TODO（rsc）：确保Ldexp可以处理边界k 
	return Ldexp(y, k)
}
